R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

source("tianfengRwrappers.R")

# library(future) 
# plan("multiprocess",workers = 8)

CA_dataset2 <- readRDS("CA_dataset2.rds")
CA_dataset1 <- readRDS("CA_dataset1.rds")
human_coronary <- readRDS("human_coronary.rds")
Idents(human_coronary) <- human_coronary$samples
human_coronary <- RenameIdents(human_coronary,'1' = 'sample1','2' = 'sample2','3' = 'sample3','4' = 'sample4')
human_coronary$samples <- Idents(human_coronary)
Idents(human_coronary) <- human_coronary$Classification1
ds2 <- readRDS("ds2.rds")

#sample info

ggsave("dataset2_sampleinfo.svg",plot = umapplot(CA_dataset2, split.by = "sample"), 
       device = svg, width = 25, height = 5)
ggsave("dataset1_sampleinfo.svg",plot = umapplot(CA_dataset1, split.by = "orig.ident"),
       device = svg, width = 15, height = 5)
ggsave("dataset0_sampleinfo.svg",plot = umapplot(human_coronary, split.by = "samples"),
       device = svg, width = 20, height = 5)

ggsave("ds2.svg",plot = umapplot(ds2), device = svg, width = 6, height = 5)
ggsave("ds1.svg",plot = umapplot(ds1), device = svg, width = 6, height = 5)
ggsave("ds0.svg",plot = umapplot(ds1), device = svg, width = 6, height = 5)

附图：所有marker基因表达热图 show 表达量最高的top5

logfc.threshold = 0.5, min.diff.pct = 0.3, pct.1 > 0.7

dataset2

CA_dataset2_markers <- FindAllMarkers(CA_dataset2, logfc.threshold = 0.5, min.diff.pct = 0.3, only.pos = T)
CA_dataset2_markers <- CA_dataset2_markers[CA_dataset2_markers$pct.1>0.7,] %>% group_by(cluster) 

genes_to_show <- CA_dataset2_markers %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_logFC)

svg(paste0("CA_dataset2_supp","_markers.svg"), height = 10, width = 15)
dhm2(CA_dataset2_markers$gene, CA_dataset2, genes_to_show$gene,"CA_dataset2_supp")
dev.off()

dataset1

CA_dataset1_markers <- FindAllMarkers(CA_dataset1, logfc.threshold = 0.5, min.diff.pct = 0.3, only.pos = T)
CA_dataset1_markers <- CA_dataset1_markers[CA_dataset1_markers$pct.1>0.7,] %>% group_by(cluster) 

genes_to_show <- CA_dataset1_markers %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_logFC)

svg(paste0("CA_dataset1_supp","_markers.svg"), height = 10, width = 15)
dhm2(CA_dataset1_markers$gene, CA_dataset1, genes_to_show$gene,"CA_dataset1_supp")
dev.off()

dataset0

human_coronary_markers <- FindAllMarkers(human_coronary, logfc.threshold = 0.5, min.diff.pct = 0.3, only.pos = T)
human_coronary_markers <- human_coronary_markers[human_coronary_markers$pct.1>0.7,] %>% group_by(cluster) 

genes_to_show <- human_coronary_markers %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_logFC)

svg(paste0("human_coronary_supp","_markers.svg"), height = 10, width = 15)
dhm2(human_coronary_markers$gene, human_coronary, genes_to_show$gene,"human_coronary_supp")
dev.off()

样本细胞比例

dataset0 冠状动脉

Idents(human_coronary) <- human_coronary$conditions
sp1 <- subset(human_coronary, idents = "sample1")
sp2 <- subset(human_coronary, idents = "sample2")
sp3 <- subset(human_coronary, idents = "sample3")
sp4 <- subset(human_coronary, idents = "sample4")
prop_mat <- cbind(prop.table(table(sp1$Classification1)),prop.table(table(sp2$Classification1)))
prop_mat2 <- cbind(prop.table(table(sp3$Classification1)),prop.table(table(sp4$Classification1)))
prop_mat <- cbind(prop_mat, prop_mat2)
colnames(prop_mat) <- levels(Idents(human_coronary))


plot_data = melt(prop_mat)
colnames(plot_data) = c('cell type','position','proportion')#修改每一列的名称

ggplot(plot_data, aes(x = `cell type`, y = proportion, fill = position)) + 
  geom_bar(stat = 'identity', position = "dodge", width=0.5) + theme_bw()

prop_plot <- ggplot(plot_data, aes(x = `cell type`, y = proportion, fill = position)) + 
  geom_bar(stat = 'identity', position = "dodge", width=0.7) + coord_cartesian(ylim = c(0,0.3))+
  theme_bw() + scale_y_continuous(expand = c(0,0)) + scale_fill_manual(values = colors_list[3:6]) +theme(
    axis.title.x = element_text(size = 15), axis.text.x = element_text(size = 15, colour = "black"),
    axis.title.y = element_text(size = 15), axis.text.y = element_text(size = 15, colour = "black"),
    legend.text = element_text(size = 20), legend.title = element_blank(), panel.grid = element_blank())

ggsave("human_coronary_prop.svg", device = svg, plot = prop_plot, width = 12, height = 6)

dataset1 颈动脉

Idents(CA_dataset1) <- CA_dataset1$orig.ident
Idents(CA_dataset1) <- c("sample1","sample2","sample3")

sp1 <- subset(CA_dataset1, idents = "sample1")
sp2 <- subset(CA_dataset1, idents = "sample2")
sp3 <- subset(CA_dataset1, idents = "sample3")
prop_mat <- cbind(prop.table(table(sp1$Classification1)),prop.table(table(sp2$Classification1)),prop.table(table(sp3$Classification1)))

colnames(prop_mat) <- levels(Idents(CA_dataset1))

plot_data = melt(prop_mat)
colnames(plot_data) = c('cell type','position','proportion')#修改每一列的名称

prop_plot <- ggplot(plot_data, aes(x = `cell type`, y = proportion, fill = position)) + 
  geom_bar(stat = 'identity', position = "dodge", width=0.7) + coord_cartesian(ylim = c(0,0.6))+
  theme_bw() + scale_y_continuous(expand = c(0,0)) + scale_fill_manual(values = colors_list[3:6]) +theme(
    axis.title.x = element_text(size = 15), axis.text.x = element_text(size = 15, colour = "black"),
    axis.title.y = element_text(size = 15), axis.text.y = element_text(size = 15, colour = "black"),
    legend.text = element_text(size = 20), legend.title = element_blank(), panel.grid = element_blank())

ggsave("CA_dataset1_prop.svg", device = svg, plot = prop_plot, width = 12, height = 6)

dataset2 颈动脉

Idents(CA_dataset2) <- factor(CA_dataset2$sample,levels = c("AC_1","AC_2","AC_3","PA_1","PA_2","PA_3"))
sp1 <- subset(CA_dataset2, idents = "AC_1")
sp2 <- subset(CA_dataset2, idents = "AC_2")
sp3 <- subset(CA_dataset2, idents = "AC_3")
sp4 <- subset(CA_dataset2, idents = "PA_1")
sp5 <- subset(CA_dataset2, idents = "PA_2")
sp6 <- subset(CA_dataset2, idents = "PA_3")

prop_mat <- cbind(prop.table(table(sp1$Classification1)),prop.table(table(sp2$Classification1)),
                  prop.table(table(sp3$Classification1)),prop.table(table(sp4$Classification1)),
                  prop.table(table(sp5$Classification1)),prop.table(table(sp6$Classification1)))

colnames(prop_mat) <- levels(Idents(CA_dataset2))

plot_data = melt(prop_mat)
colnames(plot_data) = c('cell type','position','proportion')#修改每一列的名称

prop_plot <- ggplot(plot_data, aes(x = `cell type`, y = proportion, fill = position)) + 
  geom_bar(stat = 'identity', position = "dodge", width=0.7) + coord_cartesian(ylim = c(0,0.6))+
  theme_bw() + scale_y_continuous(expand = c(0,0)) + scale_fill_manual(values = colors_list[3:8]) +theme(
    axis.title.x = element_text(size = 15), axis.text.x = element_text(size = 15, colour = "black"),
    axis.title.y = element_text(size = 15), axis.text.y = element_text(size = 15, colour = "black"),
    legend.text = element_text(size = 20), legend.title = element_blank(), panel.grid = element_blank())

ggsave("CA_dataset2_prop.svg", device = svg, plot = prop_plot, width = 16, height = 6)

XGBoost feature plot

pretrain AC–PA

fea <- read.csv("./datatable/AC_features.csv")
ggobj <- multi_featureplot(fea$Feature[1:16],ds2_AC,labels = "",label = F)
ggsave("ACpretrain_features.png", device = png, plot = ggobj, width = 8, height = 8)

fea <- read.csv("./datatable/PA_features.csv")
ggobj <- multi_featureplot(fea$Feature[1:16],ds2_PA,labels = "",label = F)
ggsave("PApretrain_features.png", device = png, plot = ggobj, width = 8, height = 8)

model AC–PA

fea <- read.csv("./datatable/ACtrain_features.csv")
ggobj <- multi_featureplot(fea$Feature[1:16],ds2_AC,labels = "",label = F)
ggsave("./supp/ACmodel_features.png", device = png, plot = ggobj, width = 10, height = 8)

fea <- read.csv("./datatable/PAtrain_features.csv")
ggobj <- multi_featureplot(fea$Feature[1:16],ds2_PA,labels = "",label = F)
ggsave("./supp/PAmodel_features.png", device = png, plot = ggobj, width = 10, height = 8)

model ds2

fea <- read.csv("./datatable/ds2_features.csv")
ggobj <- multi_featureplot(fea$Feature[1:16],ds2,labels = "",label = F)
ggsave("ds2model_features.png", device = png, plot = ggobj, width = 8, height = 8)

fea <- read.csv("./datatable/ds0_features.csv")
ggobj <- multi_featureplot(fea$Feature[1:16],ds0,labels = "",label = F)
ggsave("ds0model_features.png", device = png, plot = ggobj, width = 8, height = 8)

SMC/marker features in SMC2

Idents(ds2) <- ds2$Classification1
ds2_SMC2 <- subset(ds2, ident = "SMC2")


data2 <- FetchData(object = ds2_SMC2, vars = c("ACTA2", "TAGLN"))
rownames(data2) <-  NULL
data2$group <- "unsup"

ggplot(data2, aes(x=ACTA2, y=TAGLN, color = group, group = group)) +
  geom_point(size = 3,alpha=0.1) + 
  geom_smooth(method=lm , color="red", fill="#69b3a2", formula = 'y~x', se=TRUE) +
  theme_classic() +  theme(axis.title = element_text(size = 20,color = "black"),
        axis.text = element_text(size = 20,color = "black"),
        axis.line = element_line(size = 1),
        axis.ticks = element_line(size = 1),
        title = element_text(size = 20)) # 781 ACTA+ TAGLN+ in 792 ACTA+ 

data2 <- FetchData(object = ds2_SMC2, vars = c("SOST", "DLX5"))
rownames(data2) <-  NULL
data2$group <- "unsup"

ggplot(data2, aes(x=SOST, y=DLX5, color = group, group = group)) +
  geom_point(size = 3,alpha=0.1) + 
  geom_smooth(method=lm , color="red", fill="#69b3a2", formula = 'y~x', se=TRUE) +
  theme_classic() +  theme(axis.title = element_text(size = 20,color = "black"),
        axis.text = element_text(size = 20,color = "black"),
        axis.line = element_line(size = 1),
        axis.ticks = element_line(size = 1),
        title = element_text(size = 20))

data2 <- FetchData(object = ds2_SMC2, vars = c("DLX6-AS1", "DLX5"))
rownames(data2) <-  NULL
data2$group <- "unsup"

ggplot(data2, aes(x=`DLX6-AS1`, y=DLX5, color = group, group = group)) +
  geom_point(size = 3,alpha=0.1) + 
  geom_smooth(method=lm , color="red", fill="#69b3a2", formula = 'y~x', se=TRUE) +
  theme_classic() +  theme(axis.title = element_text(size = 20,color = "black"),
        axis.text = element_text(size = 20,color = "black"),
        axis.line = element_line(size = 1),
        axis.ticks = element_line(size = 1),
        title = element_text(size = 20))

# dim(subset(ds2_SMC2, SOST>0))[2]

dim(subset(ds2_SMC2, `DLX5`>1))[2]

dim(subset(ds2_SMC2, `DLX5`>1&SOST>1))[2]

print("...")
dim(subset(ds2_SMC2, DLX5>1))[2]

dim(subset(ds2_SMC2, `PRDM6`>1))[2]

dim(subset(ds2_SMC2, `PRDM6`>1&DLX5>1))[2]

dim(subset(ds2_SMC2, ACTA2>1))[2]

dim(subset(ds2_SMC2, TAGLN>1))[2]

dim(subset(ds2_SMC2, TAGLN>1&ACTA2>1))[2]

GO for neural progenitor in ds2

library(org.Hs.eg.db)
enrich.go <- CA_dataset2_markers[CA_dataset2_markers$cluster == "Neural progenitor",]$gene %>% enrichGO(
        OrgDb = org.Hs.eg.db,
        keyType = "SYMBOL",
        ont = "ALL", # 可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
        pAdjustMethod = "fdr",
        pvalueCutoff = 0.05,
        qvalueCutoff = 0.2,
    )
plot <- dotplot(enrich.go, title = paste("Neural progenitor", "GO"), showCategory = 15) + 
  theme_classic() + theme(text = element_text(colour = "black", size = 16), 
                          plot.title = element_text(size = 16,color="black",hjust = 0.5),
                          axis.title = element_text(size = 16,color ="black"), 
                          axis.text = element_text(size= 16,color = "black"))
ggsave("GO_Neural progenitor.svg",device = svg,height = 6,width = 10, plot = plot)

BMP receptors in different dataset

Dotplot(c("BMPR1B","BMPR1A","BMPR2","ACVR2A"),ds0)
Dotplot(c("BMPR1B","BMPR1A","BMPR2","ACVR2A"),ds1)
Dotplot(c("BMPR1B","BMPR1A","BMPR2","ACVR2A"),ds2)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.